
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The impact of divorce laws on the equilibrium in the marriage market.
% Ana Reynoso
% April 2024
%
% This file produces a panel of simulated married couples who begin their
% life under mutual consent divorce and face the introduction of unilateral
% divorce at any period after the first one.
%  
% Data: PSID 1968-1992
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%------------------------ Replication Path -------------------------------%

clc
clear all

%--- Indicate location of Replication folder:
%replication_location = 'C:\update_with_your_path';

%------------------------ Preliminaries ----------------------------------%

%--- Working directory
cd ([replication_location,'\Validation\Impact_already_married']);

%--- Path for estimates
estimates_dir = [replication_location,'\Outputs\MCD\'];

%--- Path for data
%addpath(genpath([replication_location,'\Data\Inputs']));

%--- Path for functions
addpath(genpath([replication_location,'\Internal_estimation\Inputs']));
addpath(genpath([replication_location,'\Impact_UD\Inputs']));

%--- Load estimates
load([estimates_dir, '\outputs_mcd']);
params = X;

%------------------------ Solve model ------------------------------------%
%- Behavioral parameters
parest(1,:)=params(1:9); 
parest(2,:)=params(10:18); 
sigma_f=[params(19) params(19) params(19) params(20) params(20) params(20) params(21) params(21) params(21)]; 
sigma_m=[params(22) params(23) params(24) params(22) params(23) params(24) params(22) params(23) params(24)]; 
match_qual(:,2)=params(25:33)'; 
%- Wage rate
a0m_1=params(34);
a0m_2=params(35);
a0m_3=params(36);
a0f_1=params(37);
a0f_2=params(38);
a0f_3=params(39);
%- Returns to experience
a1m_1=params(40);
a1m_2=params(41);
a1m_3=params(42);
a1f_1=params(43) ;
a1f_2=params(44) ;
a1f_3=params(45) ;
%- Returns to experience^2 and beyond
a2m_1=params(46);
a3m_1=0;
a2m_2=params(47);
a3m_2=0;
a2m_3=params(48);
a3m_3=0;
a2f_1=params(49);
a3f_1=0;
a2f_2=params(50);
a3f_2=0;
a2f_3=params(51);
a3f_3=0;
%- Returns to sahw
bm_1=params(52);
b2m_1=0;
bm_2=params(53);
b2m_2=0;
bm_3=params(54);
b2m_3=0;
bf_1=0;
b2f_1=0;
bf_2=0;
b2f_2=0;
bf_3=0;
b2f_3=0;

det_earn_m_1=[a0m_1, a1m_1, a2m_1, a3m_1, bm_1, b2m_1];
det_earn_m_2=[a0m_2, a1m_2, a2m_2, a3m_2, bm_2, b2m_2];
det_earn_m_3=[a0m_3, a1m_3, a2m_3, a3m_3, bm_3, b2m_3];
det_earn_m=[det_earn_m_1;det_earn_m_2;det_earn_m_3;det_earn_m_1;det_earn_m_2;det_earn_m_3;det_earn_m_1;det_earn_m_2;det_earn_m_3];

det_earn_f_1=[a0f_1, a1f_1, a2f_1, a3f_1, bf_1, b2f_1];
det_earn_f_2=[a0f_2, a1f_2, a2f_2, a3f_2, bf_2, b2f_2];
det_earn_f_3=[a0f_3, a1f_3, a2f_3, a3f_3, bf_3, b2f_3];
det_earn_f=[det_earn_f_1;det_earn_f_1;det_earn_f_1;det_earn_f_2;det_earn_f_2;det_earn_f_2;det_earn_f_3;det_earn_f_3;det_earn_f_3];

%- Weighted Pareto Weights
parest(3,:) = params(55:63)*0.23 + params(64:72)*0.23 + params(73:81)*0.36 +  params(82:90)*0.18;                                        

%- Decisions
[dec_mc,~, ~, ~, ~, ~, ~] = solution_backwards(n_g, T, lt, dt, H_max, disc, det_earn_f, det_earn_m, stoch_earn_f, stoch_earn_m, m_eta_f, m_eta_m, n_eta_f, n_eta_m, parest(1,:), home_prod, parest(2,:), match_qual, m_th, n_th, gamma, L_max, parest(3,:), sigma_f, sigma_m);
[dec_ud, lambda_grid, pareto_weight_index, ~, ~, ~, ~, ~, ~] = solution_backwards_ud(n_g, T, lt, dt, H_max, disc, det_earn_f, det_earn_m, stoch_earn_f, stoch_earn_m, m_eta_f, m_eta_m, n_eta_f, n_eta_m, parest(1,:), home_prod, parest(2,:), match_qual, m_th, n_th, gamma, L_max, parest(3,:), sigma_f, sigma_m);

%------------------------ Simulated housholds ----------------------------%
I = 10000;
U = 10; 
hw = ones(n_g,I,T+1,U);
lc = zeros(n_g,I,T+1,U); 
lcd=zeros(n_g,I,T-lt+1,U); 
hwd=zeros(n_g,I,T-lt+1,U); 
rp=zeros(n_g,I,T+1,U); 
for g=1:n_g
   [~, rp(g,:,:,:)]=min(abs(lambda_grid-parest(3,g)));
end
pareto_weight=NaN(n_g,I,T-lt+1,U);
pareto_weight_evol=ones(n_g,T);

seed=1812;
rng(seed); 

wft_rnd_all=rand(3,I,T); 
wft_rnd_all_rep=[wft_rnd_all(1,:,:) ; wft_rnd_all(1,:,:) ;wft_rnd_all(1,:,:) ;wft_rnd_all(2,:,:) ; wft_rnd_all(2,:,:) ;wft_rnd_all(2,:,:); wft_rnd_all(3,:,:) ; wft_rnd_all(3,:,:) ;wft_rnd_all(3,:,:) ];
wmt_rnd_all=rand(3,I,T); 
wmt_rnd_all_rep=[wmt_rnd_all(1,:,:) ; wmt_rnd_all(2,:,:) ;wmt_rnd_all(3,:,:) ;wmt_rnd_all(1,:,:) ; wmt_rnd_all(2,:,:) ;wmt_rnd_all(3,:,:) ;wmt_rnd_all(1,:,:) ; wmt_rnd_all(2,:,:) ;wmt_rnd_all(3,:,:) ];
tht_rnd_all=rand(9,I,T); 
wft_rnd=zeros(I,T); 
wmt_rnd=zeros(I,T); 
tht_rnd=zeros(I,T);
wft_gridpt=zeros(I,T);
wmt_gridpt=zeros(I,T);
th_gridpt=zeros(I,T);

for g=1:n_g
 
%-Female wages
[~,~,cdf_fu,cdf_fc]=wageproc_fn(stoch_earn_f(g,:), m_eta_f, n_eta_f);
%-Male wages
[~,~,cdf_mu,cdf_mc]=wageproc_fn(stoch_earn_m(g,:), m_eta_m, n_eta_m);
%-Theta process
[~,~,cdf_thu,cdf_thc]=thetaproc(match_qual(g,:), m_th, n_th);
%- Wage and theta nodes
wft_rnd(:,:)=wft_rnd_all_rep(g,:,:); 
wmt_rnd(:,:)=wmt_rnd_all_rep(g,:,:); 
tht_rnd(:,:)=tht_rnd_all(g,:,:);

for i=1:I   
       for t=lt:T 

           if t==lt 
               
        for f=1:n_eta_f 
            if wft_rnd(i,t)<=(cdf_fu(1,f))
               wft_gridpt(i,t)=f;
               pastnode=f;
               break
            end
        end
 
           else 
               
        for f=1:n_eta_f
            if wft_rnd(i,t)<= (cdf_fc(pastnode,f))
            wft_gridpt(i,t)=f;
            newnode=f;
            break
            end
        end
 
        pastnode=newnode;
                   
           end
           
       end
end

for i=1:I   
       for t=lt:T 
            
           if t==lt
               
        for m=1:n_eta_m 
            if wmt_rnd(i,t)<=(cdf_mu(1,m))
               wmt_gridpt(i,t)=m;
               pastnode=m;
               break
            end
        end
        
           else
               
        for m=1:n_eta_m 
            if wmt_rnd(i,t)<= (cdf_mc(pastnode,m))
            wmt_gridpt(i,t)=m;
            newnode=m;
            break
            end
        end
  
        pastnode=newnode; 
        
           end
        end
end

for i=1:I   
       for t=lt:T 
            
           if t==lt
        for th=1:n_th 
            if tht_rnd(i,t)<=(cdf_thu(1,th))
               th_gridpt(i,t)=th;
               pastnode=th;
               break
            end
        end
        
           else
               
        for th=1:n_th  
            if tht_rnd(i,t)<= (cdf_thc(pastnode,th))
            th_gridpt(i,t)=th;
            newnode=th;
            break
            end
        end
        
        pastnode=newnode; 
        
           end
           
       end
end

%- Choices
hw(g,:,lt+1:T+1,:) = 0; 
for udt=1
    for i=1:I
       for t=lt 
        lc(g,i,t,udt) = dec_mc(g,hw(g,i,t,udt), wft_gridpt(i,t), wmt_gridpt(i,t), th_gridpt(i,t), t);
        if lc(g,i,t,udt) == 2
            hw(g,i,t+1,udt) = hw(g,i,t,udt) + 1; 
        elseif lc(g,i,t,udt) == 1
            hw(g,i,t+1,udt) = hw(g,i,t,udt);
        end
       end
    
    for t=lt+1:T
        if lc(g,i,t-1,udt)==3 
            lc(g,i,t,udt)=-1;            
        elseif lc(g,i,t-1,udt)==-1 
            lc(g,i,t,udt)=-1;
            
        else 
        lc(g,i,t,udt) = dec_mc(g, hw(g,i,t,udt), wft_gridpt(i,t), wmt_gridpt(i,t), th_gridpt(i,t), t); 
        end                                     
        if lc(g,i,t,udt) == 2
            hw(g,i,t+1,udt) = hw(g,i,t,udt) + 1; 
        elseif lc(g,i,t,udt) == 1
            hw(g,i,t+1,udt) = hw(g,i,t,udt);
        else
            hw(g,i,t+1,udt) = 0; 
        end
        
    end
   end

lcd(g,:,:,udt)=lc(g,:,lt:T,udt); 
hwd(g,:,:,udt)=hw(g,:,lt:T,udt);
    
end

for udt=2:10
for i=1:I
    
    for t=lt 
        lc(g,i,t,udt) = dec_mc(g,hw(g,i,t,udt), wft_gridpt(i,t), wmt_gridpt(i,t), th_gridpt(i,t), t);
        if lc(g,i,t,udt) == 2
            hw(g,i,t+1,udt) = hw(g,i,t,udt) + 1; 
        elseif lc(g,i,t,udt) == 1
            hw(g,i,t+1,udt) = hw(g,i,t,udt);
        end
        
        rp(g,i,t+1,udt)= rp(g,i,t,udt); 
        pareto_weight(g,i,t,udt)=lambda_grid(rp(g,i,t+1,udt));
    end
   
      
     for t=lt+1:udt-1 
        if lc(g,i,t-1,udt)==3 
            lc(g,i,t,udt)=-1;            
        elseif lc(g,i,t-1,udt)==-1 
            lc(g,i,t,udt)=-1;
            
        else 
        lc(g,i,t,udt) = dec_mc(g, hw(g,i,t,udt), wft_gridpt(i,t), wmt_gridpt(i,t), th_gridpt(i,t), t); 
        end                                                                      
        if lc(g,i,t,udt) == 2
            hw(g,i,t+1,udt) = hw(g,i,t,udt) + 1; 
        elseif lc(g,i,t,udt) == 1
            hw(g,i,t+1,udt) = hw(g,i,t,udt);
        else
            hw(g,i,t+1,udt) = 0; 
        end
        
        rp(g,i,t+1,udt)= rp(g,i,t,udt); 
        pareto_weight(g,i,t,udt)=lambda_grid(rp(g,i,t+1,udt));
    end
   
    
    for t=udt:T 
        if lc(g,i,t-1,udt)==3 
            lc(g,i,t,udt)=-1;            
        elseif lc(g,i,t-1,udt)==-1 
            lc(g,i,t,udt)=-1;
            
        else 
        lc(g,i,t,udt) = dec_ud(g, rp(g,i,t,udt), hw(g,i,t,udt), wft_gridpt(i,t), wmt_gridpt(i,t), th_gridpt(i,t), t); 
        end                                                                     
        if lc(g,i,t,udt) == 2
            hw(g,i,t+1,udt) = hw(g,i,t,udt) + 1; 
            rp(g,i,t+1,udt)= pareto_weight_index(g,rp(g,i,t,udt),hw(g,i,t,udt),wft_gridpt(i,t), wmt_gridpt(i,t), th_gridpt(i,t),t); 
                                                                                                                         
            pareto_weight(g,i,t,udt)=lambda_grid(rp(g,i,t+1,udt));                                                                                                        
        elseif lc(g,i,t,udt) == 1
            hw(g,i,t+1,udt) = hw(g,i,t,udt);
            rp(g,i,t+1,udt)= pareto_weight_index(g,rp(g,i,t,udt),hw(g,i,t,udt),wft_gridpt(i,t), wmt_gridpt(i,t), th_gridpt(i,t),t);
            pareto_weight(g,i,t,udt)=lambda_grid(rp(g,i,t+1,udt));  
        else
            hw(g,i,t+1,udt) = 0; 
            rp(g,i,t+1,udt)=NaN; 
        end
        
        
    end
end 

lcd(g,:,:,udt)=lc(g,:,lt:T,udt);
hwd(g,:,:,udt)=hw(g,:,lt:T,udt);

end 

end 

%------------------------ Save outputs------------------------------------%
%--- Working directory
cd ([replication_location,'\Validation\Impact_already_married\Inputs']);

decg1=[squeeze(lcd(1,:,:,1));squeeze(lcd(1,:,:,2));squeeze(lcd(1,:,:,3));squeeze(lcd(1,:,:,4));squeeze(lcd(1,:,:,5));squeeze(lcd(1,:,:,6));squeeze(lcd(1,:,:,7));squeeze(lcd(1,:,:,8));squeeze(lcd(1,:,:,9));squeeze(lcd(1,:,:,10))];
dlmwrite('decg1.cvs',decg1);
decg2=[squeeze(lcd(2,:,:,1));squeeze(lcd(2,:,:,2));squeeze(lcd(2,:,:,3));squeeze(lcd(2,:,:,4));squeeze(lcd(2,:,:,5));squeeze(lcd(2,:,:,6));squeeze(lcd(2,:,:,7));squeeze(lcd(2,:,:,8));squeeze(lcd(2,:,:,9));squeeze(lcd(2,:,:,10))];
dlmwrite('decg2.cvs',decg2);
decg3=[squeeze(lcd(3,:,:,1));squeeze(lcd(3,:,:,2));squeeze(lcd(3,:,:,3));squeeze(lcd(3,:,:,4));squeeze(lcd(3,:,:,5));squeeze(lcd(3,:,:,6));squeeze(lcd(3,:,:,7));squeeze(lcd(3,:,:,8));squeeze(lcd(3,:,:,9));squeeze(lcd(3,:,:,10))];
dlmwrite('decg3.cvs',decg3);
decg4=[squeeze(lcd(4,:,:,1));squeeze(lcd(4,:,:,2));squeeze(lcd(4,:,:,3));squeeze(lcd(4,:,:,4));squeeze(lcd(4,:,:,5));squeeze(lcd(4,:,:,6));squeeze(lcd(4,:,:,7));squeeze(lcd(4,:,:,8));squeeze(lcd(4,:,:,9));squeeze(lcd(4,:,:,10))];
dlmwrite('decg4.cvs',decg4);
decg5=[squeeze(lcd(5,:,:,1));squeeze(lcd(5,:,:,2));squeeze(lcd(5,:,:,3));squeeze(lcd(5,:,:,4));squeeze(lcd(5,:,:,5));squeeze(lcd(5,:,:,6));squeeze(lcd(5,:,:,7));squeeze(lcd(5,:,:,8));squeeze(lcd(5,:,:,9));squeeze(lcd(5,:,:,10))];
dlmwrite('decg5.cvs',decg5);
decg6=[squeeze(lcd(6,:,:,1));squeeze(lcd(6,:,:,2));squeeze(lcd(6,:,:,3));squeeze(lcd(6,:,:,4));squeeze(lcd(6,:,:,5));squeeze(lcd(6,:,:,6));squeeze(lcd(6,:,:,7));squeeze(lcd(6,:,:,8));squeeze(lcd(6,:,:,9));squeeze(lcd(6,:,:,10))];
dlmwrite('decg6.cvs',decg6);
decg7=[squeeze(lcd(7,:,:,1));squeeze(lcd(7,:,:,2));squeeze(lcd(7,:,:,3));squeeze(lcd(7,:,:,4));squeeze(lcd(7,:,:,5));squeeze(lcd(7,:,:,6));squeeze(lcd(7,:,:,7));squeeze(lcd(7,:,:,8));squeeze(lcd(7,:,:,9));squeeze(lcd(7,:,:,10))];
dlmwrite('decg7.cvs',decg7);
decg8=[squeeze(lcd(8,:,:,1));squeeze(lcd(8,:,:,2));squeeze(lcd(8,:,:,3));squeeze(lcd(8,:,:,4));squeeze(lcd(8,:,:,5));squeeze(lcd(8,:,:,6));squeeze(lcd(8,:,:,7));squeeze(lcd(8,:,:,8));squeeze(lcd(8,:,:,9));squeeze(lcd(8,:,:,10))];
dlmwrite('decg8.cvs',decg8);
decg9=[squeeze(lcd(9,:,:,1));squeeze(lcd(9,:,:,2));squeeze(lcd(9,:,:,3));squeeze(lcd(9,:,:,4));squeeze(lcd(9,:,:,5));squeeze(lcd(9,:,:,6));squeeze(lcd(9,:,:,7));squeeze(lcd(9,:,:,8));squeeze(lcd(9,:,:,9));squeeze(lcd(9,:,:,10))];
dlmwrite('decg9.cvs',decg9);

